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Abstract 

Background: Insulin secreted by pancreatic islet /J-cells is the principal regulating hormone of glucose metabolism 
and plays a key role in controlling glucose level in blood. Impairment of the pancreatic islet function may cause 
glucose to accumulate in blood, and result in diabetes mellitus. Recent studies have shown that mitochondrial 
dysfunction has a strong negative effect on insulin secretion. 

Methods: In order to study the cause of dysfunction of pancreatic islets, a multiple cell model containing healthy 
and unhealthy cells is proposed based on an existing single cell model. A parameter that represents the function 
of mitochondria is modified for unhealthy cells. A 3-D hexagonal lattice structure is used to model the spatial 
differences among /J-cells in a pancreatic islet. The /J-cells in the model are connected through direct electrical 
connections between neighboring /J-cells. 

Results: The simulation results show that the low ratio of total mitochondrial volume over cytoplasm volume per 
/J-cell is a main reason that causes some mitochondria to lose their function. The results also show that the overall 
insulin secretion will be seriously disrupted when more than 15% of the /J-cells in pancreatic islets become 
unhealthy. 

Conclusion: Analysis of the model shows that the insulin secretion can be reinstated by increasing the 
glucokinase level. This new discovery sheds light on antidiabetic medication. 



Background 

With 25-30 percent of adults in the developed world hav- 
ing a high risk of diabetes, and 24 million people in the 
United States with diabetes, diabetes mellitus ranks as a 
principal cause of death [1]. The incidence of diabetes 
mellitus is expected to double within the coming 20 years. 
After a meal, healthy individuals secrete insulin into the 
bloodstream signaling the consumption of glucose pro- 
duced from the digested food. Thus insulin has an impor- 
tant signaling role in the maintenance of low blood 
glucose levels via glucose consumption. Type II diabetes 
may result from disruption of this signaling process [2,3]. 
Hence understanding the insulin secretion mechanism is 
vitally important. 
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The main pathways involved in glucose-stimulated insu- 
lin secretion are shown in Figure 1. Glucose is imported 
into /J-cells and undergoes the glycolysis process. Three pri- 
mary products, pyruvate, malate, and nicotinamide adenine 
dinucleotide plus hydrogen (NADH), are produced in gly- 
colysis [4] and transported into mitochondria. In mitochon- 
dria pyruvate and malate are consumed in the tricarboxylic 
acid (TCA) cycle to produce fumarate. Adenosine tripho- 
sphate (ATP), the most important energy resource in the 
cell, is produced from fumarate and NADH through oxida- 
tive phosphorylation and released to the cytoplasm [5]. The 
increasing ATP level in the cytoplasm blocks the ATP-sen- 
sitive K + channel [6-8], resulting in a higher intracellular IC 
level and depolarization of the cell membrane. This results 
in the opening of voltage-dependent Ca + channels. Cyto- 
plasmic Ca + is transported into the /J-cell from the extra- 
cellular fluid through the voltage-gated Ca + channels 
[9-11]. When both the cytoplasmic ATP and Ca 2+ levels 
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rise high enough, the vesicles containing insulin fuse with 
the cell plasma membrane, releasing insulin from the cell. 
The increasing formation and export of ATP to the entire 
cell will raise the Ca 2+ levels in the mitochondria as well. 
High mitochondrial Ca 2+ levels eventually collapse the 
mitochondrial membrane potential, shutting off ATP pro- 
duction, possibly by opening a low-conductance state of 
the permeability transition pore (PTP) in the mitochondrial 
membrane [12-14]. The rise and fall in the mitochondrial 
membrane potential and intramitochondrial Ca + concen- 
tration is the mitochondrial oscillation. The frequency of 
this oscillation depends on how rapidly the mitochondria 
can build up its membrane potential and transport Ca 2+ . 
Besides the metabolic pathway in a single /J-cell, cells in 
each pancreatic islet are also connected through direct elec- 
trical connections between neighboring /J-cells. This con- 
nection is sufficient to synchronize both electrical bursting 
activity and metabolic oscillations [15]. 

The cells, in fact, have differing capabilities in consum- 
ing glucose, since their differing physical properties cause 
different rates of electron transport chain function, ATP 
production, Ca 2+ uptake and release, and cell membrane 
potential. Furthermore, some /J-cells do not function nor- 
mally as just described. The term unhealthy is used to 
describe cells that cannot consume glucose normally. 
Unhealthy cells interact strongly with adjacent healthy 
cells, because of their electrical connections. The total 
insulin secretion is the aggregation of the insulin secreted 
by all pancreatic islet cells. Its dynamics may be destroyed 
by a small fraction of unhealthy cells even though a large 
majority of the cells are still healthy. Consequently, under- 
standing how a cohort of unhealthy cells can affect the 
insulin secretion process in other healthy cells is very 
important, and can lead to a better understanding of the 
cause of diabetes. 

This present work proposes a mathematical model to 
study scenarios with unhealthy cells leading to insulin 



oscillation failure. The new model is based on a mathe- 
matical insulin secretion model of the oscillation and 
metabolic pathway of insulin secretion in a single /J-cell, 
due to Bertram et al. [16-18]. Bertram's original model is 
extended to model spatial and coupling effects [14] due 
to the electrical connections between neighboring /J-cells. 
Mitochondrial dysfunction is believed to have a strong 
negative effect on insulin secretion [19-21], and also 
mitochondrial dysfunction is believed to be a main cause 
of unhealthy cells, based on the recent data of [22-24]. 
The ratio of the total mitochondrial volume to the cyto- 
plasm volume per /J-cell, a parameter corresponding to 
mitochondria function, is used to define unhealthy cells, 
as distinguished from healthy cells. The effect on the islet 
insulin oscillation of unhealthy cells coupled with healthy 
cells is then studied in detail. Multiple cell simulations 
demonstrate insulin oscillation malfunction when the 
fraction of unhealthy cells exceeds approximately 15%. 
To verify this result, a simplified coupling topology is 
used to study the effect of one unhealthy cell on neigh- 
boring healthy cells. The latter study confirmed the pre- 
vious more complicated simulation. 

In order to understand the cause of insulin secretion 
failure resulted from unhealthy cells, an eight-cell model 
is built and studied in details. We discovered a critical 
difference of the dynamics between healthy cells and 
unhealthy cells. Our simulation results demonstrate that 
stimulating glucokinase can make unhealthy pancreatic 
islets function normally. Based on the discovery, a possi- 
ble strategy for antidiabetic medicine is proposed. Our 
strategy is consistent with recent antidiabetic medicine 
development [25-30] that identifies glucokinase as a 
major drug target. 

Methods 

Single /J-cell modeling 

The mathematical model of a single pancreatic /J-cell is 
based on the deterministic model introduced by Bertram 
et al. [16-18]. The model has four components: glycolytic, 
mitochondrial metabolism, cytoplasmic intermediate, and 
plasma membrane. Each component of the model is con- 
nected by variables shared between components as illu- 
strated in Figure 2. 

The glycolytic component of the model is a simplified 
model of the glycolysis pathway [31] that converts glu- 
cose into pyruvate. This pathway generates oscillations 
via a critical enzyme phosphofructokinase (PFK), the rate 
limiting variable in glycolysis and an important control. 
PFK is allosterically inhibited by ATP [32]. As a result 
high ATP levels will cause glycolysis to slow down. PFK 
phosphorylates fructose-6-phosphate (F6P) to fructose- 
1,6-bisphosphate (FBP), which is a substrate for glycera- 
dehyde 3-P dehydrogenase (GPDH). The GPDH reaction 
is in a rapid equilibrium, therefore used as the input to 
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Figure 2 Model components and interconnections. The whole model is divided into four components: glycolytic, mitochondrial metabolism, 
cytoplasmic intermediate, and plasma membrane. 



the mitochondrial component, replacing pyruvate dehy- 
drogenase (PDH) in the model. The glycolysis model 
consists of equations for the concentrations [G6P] and 
[FBP]. Here [G6P] is the glucose 6-phosphate concentra- 
tion, which is assumed to be in rapid equilibrium with 
[F6P], and the /s are fluxes. 



d[G6P] 
dt 



'■ JGK — JPFK, 



d[FBP] 1 

-, = JpFK — —Jc.PDH- 

dt 2 



(la) 



(lb) 



The model for the mitochondrial component is based on 
the detailed Magnus-Keizer model of mitochondrial meta- 
bolism [33]. The tricarboxylic acid cycle (TCA cycle) and 
oxidative phosphorylation occur in mitochondria. The 
TCA cycle is a series of enzyme-catalyzed chemical reac- 
tions whose main role is to convert pyruvate to NADH 
and reduced flavin adenine dinucleotide (FADH), and 
then feed them into oxidative phosphorylation, which uses 
energy released by the oxidation of nutrients to produce 
ATP [34]. The model includes an expression for Ca 2+ - 
dependent dehydrogenases of the citric acid cycle, yielding 
NADH. The NADH supplies electrons for the electron 
transport chain, which produces a membrane potential Ay/ 
across the mitochondrial inner membrane. The flux of 
protons down the electrical gradient through the F1F0- 
ATP synthase converts ADP to ATP. Finally, Ca + enters 
the mitochondria through Ca + uniporters and is trans- 
ported out through Na + /Ca + exchangers [35]. The mito- 
chondrial compartment has variables for the NADH, 
ADP, and mitochondrial Ca 2+ concentrations, and for Ay/ 
The ATP concentration is calculated from the ADP con- 
centration by assuming that the sum of the two is con- 
served. Differential equations for the mitochondrial 
variables are 



d[NADH„ 



dt 



yVdh-Jo), 



(2a) 



d _ 
dt 

d[Ca m ] 
dt 

d[ADP, 
dt 



(/h,R«s — Jn.alp ~ JANT ~ JhMIi — JNaCa ~ 2Juni)/C m , (2b) 
= fm {Juni — Jn« Ca ) i (^c) 



= y{Ja 



Jfifo)- 
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There are five variables in the cytoplasmic intermediate 
and plasma membrane components, which describe the 
plasma membrane potential V and the cytosolic ADP con- 
centration, the fraction n of the opened delayed rectifying 
K + channels and cytosolic Ca 2+ concentration, as well as 
the ER Ca + concentration. The differential equations are 



dV 
dt 



{Ik + ha + h(Ca) + Ik{atp))/C, 

dn noo(V) — n 
dt x 

d[Ca c ] 



dt 

d[Ca er ] 
dt 

d[ADP c 
dt 



1 fcijmem + Jer + ^/m)/ 
= -feriVc/VerVer, 

= Jhyd — KJANT, 



(3a) 
(3b) 
(3c) 
(3d) 
(3e) 

are the ionic current on 



where I K , I Ca , I K(Ca ), and I k(AT p) 
the membrane, C is the membrane capacitance, f c is the 
fraction of free Ca 2+ in the cytosol, } mem is the flux of Ca 2+ 



Pu et al. BMC Medical Genomics 201 3, 6(Suppl 3):S6 
http://www.biomedcentral.eom/1755-8794/6/S3/S6 



Page 4 of 13 



across the plasma membrane, J m is the flux of Ca + out of 
the mitochondria, which is scaled by the mitochondria/ 
cytosol volume ratio K, J er is the flux out of the ER,^ r is the 
fraction of free Ca 2+ in the ER, V c and V er are the volumes 
of the cytosolic and ER compartments, respectively, n is the 
fraction of the open delayed rectifying K + channels, r is the 
relaxation time for the open and close reactions of the 
delayed rectifying K + channels to reach equilibrium, and 
the steady state function for n is « M (V) = (1 + exp(-(V + 
l6)/5))-\ For more details on these terms in the differential 
equations of the model, refer to Bertram et al. 

Multiple /J-cells modeling 

There are about 1, 000 /J-cells in each pancreatic islet. 
They interact with neighboring /3-cells through direct 



electrical connections. These interactions are modeled 
with a change in (3a): 



dVj 

-jj- = —(Ik + lea + Ik{Co) + Ik{atp) + Ij)/C, 

ieN(j) 



(4) 



where Ij is the current for cell /, gc is the electrical cou- 
pling conductance, and N (/') is the index set for all neigh- 
bor cells of cell /. The neighbors surrounding a cell / are 
detected by considering the position of a cell / within a 
3-D hexagonal lattice (different from the 3-D von Neu- 
mann neighborhoods of cellular automata theory). As 
shown in Figure 3, each internal cell is connected with six 



Level 1: 



Level 2: 



Level 3: 




Figure 3 Three layers of jS-cells in 3-D hexagonal lattice 
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cells in the same layer, three cells in the upper layer and 
three other cells in the lower layer. For example, cell 23 in 
Figure 3 is connected with cells 18, 19, 22, 24, 27, 28 in its 
own layer, and connected with cells 7, 11, 12 in the upper 
layer, and cells 34, 35, 39 in the lower layer. Face, edge, 
and corner cells have between three and nine neighbors, 
depending on their location. Table 1 is used to calculate 
the coordinates of the neighbors of cell {a, b, c) in the 3-D 
hexagonal lattice. 

Besides the spatial differences, each cell has its own 
parameter set. Among all the parameters of the single 
cell model, k, the ratio of mitochondria volume to cytosol 
volume, represents the function of mitochondria. The 
single cell study has shown that the oscillation behavior 
of a cell is very sensitive to k, and thus k is used to 
describe the cell differences. The value of n is different 
for each cell in the 3-D structure, while all other cell 
parameters are set the same. 

Results and discussion 

The differential equations for the model, implemented in 
Fortran 90, were solved using LSODE in ODEPACK [36]. 
Generated data were plotted using Matlab. The simulations 
are performed on a Mac OS 2.4 GHz Intel(R) Core 2 Duo 
CPU with 4GB memory. 

Single /J-cell simulation 

The simulation model for multiple cells has at its core the 
single cell model. Since the total amount of ATP and ADP 
(in the single cell model) is conserved, only ADP appears 
in the differential equation 

- -J ^ = Jhyd - KJANT, (5) 

where /w is the hydrolysis rate of ATP to ADP, k (k = 
0.0733 in the standard model [37]) is the ratio of the total 
functional mitochondrial volume to the cytoplasm volume, 



Table 1 The coordinates of the neighbors of the cell with 
coordinates (a, b, c) 



Layer 


Coordinates (x, y, z) 




(a - 1, fa, c) 


Upper Layer 


(a - 1, b + 1, c+ 1) 




{a - 1, b, c+ 1) 




{a, b, c- 1) 




(a, fa, c + 1) 


Same Layer 


(a, fa - 1, c - 1) 




(a, fa - 1, c) 




(a, fa + 1, c) 




(a, fa + 1, c + 1) 




(a + 1, b - 1, c - 1) 


Lower Layer 


(a + 1, fa, c - 1) 




(a + 1, fa, c) 



and Jant is the flux through the adenine nucleotide trans- 
locator, which exchanges ADP and ATP molecules 
between the cytoplasm and the mitochondria. With k = 
0.0733, both the cell membrane potential and insulin 
secretion show periodic bursts of rapid oscillation, as illu- 
strated in Figure 1 of Additional file 1. The slowest com- 
ponent of the compound bursting is due to oscillatory 
glycolysis, reflected by an oscillatory FBP concentration. 
The oscillatory glycolysis causes slow oscillations in ADP, 
which superimpose with the faster ADP oscillations driven 
by Ca 2+ . The multimodal ADP rhythm leads to oscillations 
in the conductance of the ATP-dependent potassium 
channel, which drives the burst episodes of the membrane 
potential V, which then results in compound bursting of 
intracellular calcium, leading to pulsatile insulin secretion. 

In truth, though, k varies between cells; specifically, cells 
with some mitochondria weakened due to aging will 
effectively have a k smaller than the default value. It is 
therefore paramount to study the sensitivity of the cell 
membrane potential V and insulin secretion / patterns to 
the parameter k. For this work, Bertram's single cell 
model is modified by varying the parameter k to represent 
different levels of functional mitochondria within different 
cells. In terms of changes in the patterns of V, Ca 2+ , and /, 
the bursting behavior of the model is classified into four 
categories as a function of the volume of functioning mito- 
chondria in the /J-cell, in order of occurrence: burst forma- 
tion, periodic burst, burst loss, and decoupling. 

The bursting pattern starts to form at around 85% of the 
regular mitochondrial volume, which is when the total 
volume of the mitochondria is about 6% of the cytoplas- 
mic volume. If the volume of functioning mitochondria is 
lower than this limit, insufficient ATP is produced to 
allow the insulin burst. In the burst formation category 
(Additional file 1: Figure 2) each burst consists of a small 
number of spikes in plasma membrane voltage, cytoplas- 
mic calcium ion concentration, and the insulin secretion 
rate. If the mitochondrial volume fraction is increased 
beyond the value used in the standard model, the behavior 
of the burst begins to change in pattern, entering the burst 
loss category (Additional file 1: Figure 3). In the burst loss 
category every other burst decreases in duration until it is 
lost completely, while the remaining bursts lengthen their 
duration. The net effect is a decrease in the fraction of 
time spent in each burst. 

When the volume fraction of the mitochondria rises 
over about 0.095 the oscillation patterns undergo a radical 
alteration (Additional file 1: Figure 4). The frequency of 
the cell membrane voltage oscillation within the bursts 
becomes very large. While the membrane voltage con- 
tinues to undergo bursts of rapid oscillations separated by 
periods of no oscillation, this no longer drives the calcium 
ion concentration, which then normally drives the insulin 
secretion. Thus this behavior will be called the decoupling 
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category. While insulin release does undergo a long and 
slow periodic variation, the rapid spikes that characterize 
the normal burst no longer occur. 

In order to quantify the development and eventual loss 
of the bursting behavior, define the "bursting fraction" as 
the portion of the total time taken up by the bursts, where 
a burst duration is defined as the time from the first to the 
last peak in a burst. Bertram's single cell original model 
has a bursting fraction of approximately 0.55 (meaning 
that a burst in the model variables' oscillations is occurring 
for 55% of the time). When the cell mitochondrial volume 
fraction is changed, this bursting fraction changes mark- 
edly (cf. Figure 4), and the patterns of the calcium ion con- 
centration Ca 2+ , the cell membrane potential V, and the 
insulin secretion rate / all alter noticeably. Figure 4 shows 
the ranges of the mitochondrial volume ratios correspond- 
ing to the four classifications: burst formation, periodic 
bursts, burst loss, and decoupling. Since the classifications 
are qualitative, the divisions between the classifications are 
blurred and the ranges are drawn overlapping. The mem- 
brane voltage, calcium ion concentration, and insulin 
release all have nearly identical bursting fractions for the 
first three categories (formation, periodic busts, and loss). 
In the decoupling category the bursting fractions in the 
calcium ion concentration and insulin release rapidly fall 



to zero, while the voltage bursting fraction holds relatively 
constant. The normal bursting behavior of the single 
/3-cell model occurs only within a narrow range of mito- 
chondrial volumes, from roughly 7% to 8% of the cellular 
volume. 

The reason for the absence of the bursting function for 
low mitochondrial volume is straightforward. With a 
lower volume of functioning mitochondria the cytoplasmic 
ATP levels are lower and are insufficient to trigger the 
remaining stages of the signaling pathway for insulin 
release (Figure 1). The change in behavior at higher mito- 
chondrial volumes, where cytoplasmic ATP levels are 
slightly higher, has a more complex mechanism. This 
model of the insulin secretion is built around a fast central 
oscillator [10,38,39] consisting of the plasma membrane 
potential V and the fraction n of open delayed rectifying 
potassium ion channels, denoted by K + in Figure 1. As the 
symbol implies, the opening and closing of these potas- 
sium channels is controlled by the cytoplasmic ATP con- 
centration. During the bursts these two variables execute a 
fast oscillation (Figure 5), which drives the other oscilla- 
tions of this model system. Between bursts, this oscillation 
stops and the V and n values move over to the "tail of the 
Q" in Figure 4, until the next burst begins. The period 
of this central oscillator driving the whole model is 
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Figure 4 Bursting fraction for the cell membrane potential V , the calcium ion concentration Ca 2+ , and the insulin secretion rate /. The 

X-axis represents the ratio of functional mitochondrial volume and the cytoplasmic volume. The Y-axis represents the bursting fraction. 
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controlled by the ATP dependence on the opening of the 
potassium channels. Higher ATP levels lead to a faster 
oscillation of V and n, which (in this model) becomes too 
fast for the dynamics of the cytoplasmic calcium ion con- 
centration to follow, decoupling the rest of the system 
from the central oscillator. Therefore, altering the mito- 
chondrial/cytosol volume ratio in the model, which corre- 
sponds to the amount of functional mitochondria in the 
cell, alters the amount of ATP production and through 
this the level and pattern of insulin secretion. In order to 
reflect cell differences, the value of k is different for each 
cell in the 3-D structure, with all other cell parameters 
being the same. 

The insulin bursts are totally absent when the volume 
fraction of the mitochondria is either less than 0.06 or 
greater than 0.095, which suggests an alternative definition 
for "healthy cell" and "unhealthy cell". Define a healthy 
cell as having a mitochondria volume fraction k between 
0.06 and 0.095, and a cell with k out of this range as an 
unhealthy cell. It is plausible that some mitochondria lose 
their function, effectively resulting in a smaller k. The 
ensuing numerical experiments use k = 0.05 for unhealthy 
cells and the standard value k = 0.0733 for healthy cells. 

Multiple /J-cells simulation 

Although there are about 1,000 cells in each pancreatic 
islet, for multiple /J-cells simulations, consider first the 
case with 125 cells coupled spatially in a 3-D hexagonal 
lattice. The justification for using 125 rather than 1000 
cells is a pragmatic one-the CPU time for simulating 
1,000 cells is rather long. In a 1,000 cell heterogeneous 
model, each single cell model has ten variables, yielding a 
10,000-dimensional ODE, which required 26 hours to 
solve till time t = 2 x 10 6 . Furthermore, the oscillation pat- 
terns observed from 125 cells are qualitatively very similar 
to those observed from 1,000 cells. 

If there are no unhealthy cells at all, the membrane 
potentials of all the cells synchronize after the coupling is 
turned on at time 400,000 milliseconds (ms) by changing 
gc from 0 to 150. This simulation of 125 cells without 
unhealthy cells is shown in Figures 1 (membrane potential) 
and 2 (total insulin secretion) in Additional file 2. The 
curves in the membrane potential plot are out of phase at 
time t = 0, but soon after 400,000 ms, these curves coalesce 
(see Additional file 2: Figure 1). Because the insulin levels 
of some cells are high while those of other cells are low, 
the total insulin is relatively flat before synchronization. 
Immediately after the coupling is turned on, the total insu- 
lin secretion shows bursts and its value rises to a hundred 
times that of a single cell, because there are more than a 
hundred cells synchronized and releasing insulin in phase. 

Focus on total insulin secretion to see how unhealthy 
cells, through the 3-D coupling in the hexagon structure, 



affect the total insulin secretion. To save computational 
time the coupling is turned on at the beginning of the 
simulations, t = 0. Figure 3 in Additional file 2 shows the 
resulting total insulin behavior with 10% of the cells 
being unhealthy spread uniformly in the 3-D hexagonal 
structure. The total insulin, as in the case of 100% healthy 
cells, shows periodic oscillations and maintains a reason- 
able level. When the percent of cells being unhealthy 
increases to 15%, the oscillations of total insulin still look 
normal (see Additional file 2: Figure 4), but now some 
bursts have fewer spikes. As the percentage of cells being 
unhealthy increases to 20% and 30% from 10% and 15% 
of cells being unhealthy, the spikes within each burst 
become much less numerous (shown in Additional file 2: 
Figures 5 and 6). These bursts are also much more irre- 
gular, and even more significantly, totally disappear after 
2.25 x 10 6 ms (Additional file 2: Figure 6). In summary, 
the cohort of unhealthy cells dominates the global beha- 
vior, resulting in a level of total insulin too low to main- 
tain proper pancreatic islet function. The conclusion is 
that if there is more than approximately 15% of cells 
being unhealthy, the function of the pancreatic islet will 
be severely affected. 

To verify the conclusions from simulations with 125 
cells, simulations were also performed for the model with 
1,000 cells. This 1000 cell model cannot be run to as long 
a time as the 125 cell model, but observe that when, for a 
certain time, insulin secretion cannot generate enough 
spikes, the pancreatic islet can be considered to be mal- 
functioning. Hence the simulation monitors the numbers 
of spikes in the last five bursts, and if the mean number of 
spikes is below three, deems that the overall system is mal- 
functioning and halts. The simulations for 1,000 cells with 
30%, 20%, 15% of cells being unhealthy, all considered 
malfunctioning systems, are shown in Figures 7, 8, and 
9 of Additional file 2 respectively. A simulation for 10% of 
cells being unhealthy found no malfunction in an extre- 
mely long time (comparable to the time for the 125 cell 
model runs). In summary, the 125 cell and 1000 cell simu- 
lations yield similar conclusions: If the percentage of cells 
being unhealthy is larger than approximately 15%, the sys- 
tem will malfunction; at 15% of cells being unhealthy, the 
system still functions but is close to malfunctioning; below 
10% of cells being unhealthy, the system can definitely 
function very well. 

Another ratio, the number of links between an 
unhealthy cell and a healthy cell to the total number of 
links between all cells, is also calculated for each simula- 
tion; the value of this ratio is quite close to the number 
of unhealthy cells to total number of cells ratio, see 
Table 2. Therefore, for further analysis, using 125 cells 
and the ratio of number of unhealthy cells to total num- 
ber of cells should suffice. 
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Simplified multiple /J-cells simulation 

In order to understand further how unhealthy cells affect 
the healthy cells, consider a simplified case with only one 
unhealthy cell in the 3-D structure, as shown in Figure 6. 
The red dot represents the unhealthy cell, the blue dots 
healthy cells. The unhealthy cell is connected with all the 
healthy cells, while the healthy cells themselves are con- 
nected through three rings corresponding to three layers 
in the 3-D structure. The experiment starts with two 
cells, one healthy and one unhealthy. Then the number 
of healthy cells is increased one by one. If the total num- 
ber of cells is smaller than 13, the corresponding sub- 
graph of Figure 6 is used. Figure 1 in Additional file 3 
shows the total insulin of one unhealthy cell and one 
healthy cell (coupling starts at 3 x 10 5 ms). According to 
the "malfunction criterion" (the average number of spikes 
in the last five bursts is smaller than three), the system 



Table 2 Comparison of Two Different Ratios 



Ratio of Links 


Ratio of Cells 


0.0945 


0.10 


0.1477 


0.15 


0.1963 


0.20 


0.2994 


0.30 



malfunctions. Figure 2 in Additional file 3 shows the 
total insulin for one unhealthy cell and three healthy 
cells. The simulations show that when the number of 
healthy cells reaches three (marginally) or four, the total 
insulin shows functioning bursts. This suggests that one 
unhealthy cell needs at least four healthy cells to make 
the whole system function well based on the topology in 
Figure 6. To study the case of two unhealthy cells with 
exactly the same topology and connected to the same 
number of healthy cells, the topology in Figure 7 is used. 
The number of healthy cells needed to rescue two 
unhealthy cells is eight (marginally) or nine, which is con- 
sistent with the result of one unhealthy cell in Figure 6. 
Therefore one unhealthy cell needs at least four healthy 
cells to make the overall system function well, i.e., the sys- 
tem probably can tolerate at most 20% of cells being 
unhealthy. This result matches with the previous simula- 
tions with 125 and 1,000 cells in the 3-D topology. 

The cause of insulin secretion failure 

To study possible reasons for the failure of insulin secre- 
tion caused by unhealthy cells, consider a 2 x 2 x 2 
model of eight cells, in which three unhealthy cells are 
uniformly distributed. Figures in Additional file 4 show 
the behaviors of all the variables in this model for each of 
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Figure 6 Topology for one unhealthy cell. A simplified case with only one unhealthy cell in the 3-D structure. The red dot (center dot) 
represents the unhealthy cell, the blue dots are healthy cells. 



the eight cells in each figure. The insulin secretion starts 
to fail at around time 2 x 10 ms as shown in Figure 8. 
Except for the plot of [G6P] (see Additional file 4: Figure 
10), it is hard to distinguish healthy cells from unhealthy 
cells in all the plots, such as that of membrane potential 
in Figure 5 of Additional file 4. In Figure 10 of Additional 
file 4 eight curves are partitioned into two groups. The 
upper group consists of healthy cells, while the lower one 
consists of unhealthy cells. These two groups of curves 
do not intersect with each other after the failure, but they 
are mixed together before the failure. This observation 
suggests that G6P might be closely related to the insulin 
secretion failure. In order to test this hypothesis, the two 
groups of curves were manually reset to see whether 
insulin secretion can be reinstated: if the insulin secretion 
fails for 5 x 10 ms, the values of [G6P] in all cells are 
reset to 279, the initial value of [G6P] used in the simula- 
tion. The oscillations of insulin secretion are resumed 
(see Additional file 5: Figures 1 and 2). The black vertical 
lines in these figures are the times when the values of 
[G6P] are equalized. The results show that the divergence 
of the two groups of [G6P] leads to the disappearance of 
insulin secretion oscillations. 

Since an unhealthy cell has a smaller mitochondria/ 
cytosol volume ratio than that of a healthy cell, unhealthy 



cells usually have lower ATP (from oxidative phosphoryla- 
tion) levels. ATP has a negative feedback to phosphofructo- 
kinase (PFK) in glycolysis. When the ATP level is lower, 
the PFK reaction becomes faster, which directly consumes 
G6P faster leading to lower G6P levels in unhealthy cells 
than in healthy cells. As shown in Figure 10 of Additional 
file 4 the levels of G6P in healthy cells are higher than 
those in unhealthy cells. At the same time, because of the 
low production rate of ATP in unhealthy cells, the ratio of 
ATP to ADP gets lower. ATP-sensitive K + channels in the 
plasma membrane are activated by ADP and inactivated by 
ATP, so the ratio of these nucleotides determines the frac- 
tion of open ATP-sensitive K + channels. When the ATP/ 
ADP ratio is low there is an increase in the number of 
open ATP-sensitive K + channels, which results in the diffi- 
culty of membrane depolarization. Voltage-dependent Ca + 
channels are blocked. Since the insulin is secreted when 
Ca 2+ exceeds a certain level, the blocked Ca 2+ channels will 
reduce insulin secretion. Therefore, raising the levels of 
G6P in unhealthy cells back to normal will bring back nor- 
mal insulin secretion. This analysis suggests the hypothesis 
that increasing the value of any substance ahead of the 
ATP synthesis in the pathway, such as the substances FBP 
and NADH in mitochondria, reinstates the oscillations. 
Manually resetting the values of either [FBP] or [NADH] 
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in all the cells to the same initial value confirms the conjec- 
ture: not only can G6P restart the insulin oscillations, but 
also FBP and NADHm can restart the oscillations. 

Since G6P and FBP are both substances in the glycolysis 
pathway, if the rate of glycolysis were faster, the oscillations 
of insulin might be resumed as well. In order to test this 
conjecture, the glucokinase level is raised fourfold after the 
insulin failure is detected. Glucokinase is the enzyme that 
phosphorylates glucose to glucose-6-phosphate (G6P). 
Figures 3 and 4 in Additional file 5 show that the oscilla- 
tions of insulin are resumed when the glucokinase level 
increased. Two oscillation periods are shown in Figure 5 of 
Additional file 5 (from the earlier 125 cell model): The 
longer period (slow oscillations) is caused by the glycolytic 
oscillations, while Ca 2+ feedback is responsible for the fast 
oscillations. Only the short period remains in Figure 6 of 
Additional file 5. This observation suggests that if the 
longer period could be reestablished in Figure 6 of Addi- 
tional file 5 the insulin oscillations might be reinstated. 
Since the glycolysis pathway is responsible for the slow 
oscillations, that pathway may need some stimulus. The 
simulation results in Figures 4 and 5 of Additional file 5 
demonstrate that by increasing the glucokinase level in the 
glycolysis pathway, the insulin oscillations can be resumed 
after the failure. 

Figure 7 in Additional file 5 shows the simulation result 
when the three unhealthy cells are placed together in the 
2x2x2 grid. With a high ratio of unhealthy cells to total 



cells (more than .20), the insulin secretion fails. If the glu- 
cokinase level in the unhealthy cells is raised fourfold and 
then kept unchanged during the simulation, the total insu- 
lin secretion becomes normal (see Additional file 5: Figure 
8). From Figure 9 in Additional file 5 all the unhealthy 
cells look like the healthy cells. Even though those 
unhealthy cells are still unhealthy, if their glucokinase 
levels are high enough, the overall behavior of the system 
will be dominated by the healthy cells. If there is no 
healthy cell, only unhealthy cells with high glucokinase 
levels, the system can't be repaired, see Figure 9. There- 
fore, the unhealthy cells are rendered harmless by high 
glucokinase level and other healthy cells in the network. 

Conclusion 

In conclusion, the insulin secretion of pancreatic islets 
usually will be functionally destroyed when there are more 
than 20% of cells being unhealthy among all cells. The 
more unhealthy cells there are, the more irregular insulin 
secretion will be. Increasing the level of glucokinase can 
make the pancreatic islet function normally when there is 
a high fraction of unhealthy cells by increasing the glucose 
absorption of the glycolysis pathway. This has implications 
for the clinical treatment of type II diabetes. Currently 
there are three classes of medications used to treat type II 
diabetes. The first treatment is to increase the amount of 
insulin secreted by the pancreas by inhibiting the opened 
delayed rectifying K + channels. The second treatment is to 




Time X 10 5 

Figure 9 Total insulin secretion of the eight unhealthy cells system with the stimulation on glucokinase of the cells. All the cells are 
placed in the 2x2x2 topology. 
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increase the sensitivity of target organs to insulin. The 
third treatment is to decrease the rate at which glucose is 
absorbed from the gastrointestinal tact, which is a method 
to reduce the glucose uptake from food. It appears that all 
the available treatments are insufficient to stem the tide. 
Therefore, new treatments are currently under investiga- 
tion including the development of therapeutic agents with 
novel action mechanisms. Recently, researchers have iden- 
tified glucokinase as an outstanding drug target for devel- 
oping antidiabetic medicines [25-30]. Assuming the 
mathematical models are valid, the simulation results 
demonstrate that stimulating glucokinase can make 
unhealthy pancreatic islets function normally, consistent 
with new antidiabetic medicines. Such multiple cell mod- 
els are good candidates for guiding the development of the 
next generation of antidiabetic medicines. 

Additional material 
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